# With Village Fixed Effects ---------------------------------------------------

# Remove observations that contain missing values for electricity supply dimensions
pool_elec <- pool %>% dplyr::select(elec_respvars, village, hh, age, hhsize, bplaay,
                                    scst, sc, st, wave, gender, religion, education, lnexp,
                                    weights) %>% 
  na.omit()


# Duration Day (Total hours per day)
model.dd.vfe.duration_day.scst = fixest::feols(duration_day ~ scst*wave + age + gender +
                                                 religion + education + hhsize +
                                                 lnexp + bplaay | village, 
                                               data = pool_elec, weights = pool_elec$weights)

# Duration Night (Hours in the evening)
model.dd.vfe.duration_night.scst = fixest::feols(duration_night ~ scst*wave + age + gender +
                                                   religion + education + hhsize +
                                                   lnexp + bplaay | village, 
                                                 data = pool_elec, weights = pool_elec$weights)

# Reliability (Outage days per month)
model.dd.vfe.reliability.scst = fixest::feols(reliability ~ scst*wave + age + gender +
                                                religion + education + hhsize +
                                                lnexp + bplaay | village, 
                                              data = pool_elec, weights = pool_elec$weights)

# Voltage Fluctuation (Voltage fluctuation days per month)
model.dd.vfe.quality_fluc.scst = fixest::feols(quality_fluc ~ scst*wave + age + gender +
                                                 religion + education + hhsize +
                                                 lnexp + bplaay | village, 
                                               data = pool_elec, weights = pool_elec$weights)

# Low Voltage (Low Voltage days per month)
model.dd.vfe.quality_low.scst = fixest::feols(quality_low ~ scst*wave + age + gender +
                                                religion + education + hhsize +
                                                lnexp + bplaay | village, 
                                              data = pool_elec, weights = pool_elec$weights)

# With Household Fixed Effects -------------------------------------------------

# Duration Day (Total hours per day)
model.dd.hfe.duration_day.scst = fixest::feols(duration_day ~ scst*wave + hhsize +
                                                 lnexp + bplaay | hh, 
                                               data = pool_elec, weights = pool_elec$weights)

# Duration Night (Hours in the evening)
model.dd.hfe.duration_night.scst = fixest::feols(duration_night ~ scst*wave + hhsize +
                                                   lnexp + bplaay | hh, 
                                                 data = pool_elec, weights = pool_elec$weights)

# Reliability (Outage days per month)
model.dd.hfe.reliability.scst = fixest::feols(reliability ~ scst*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool_elec, weights = pool_elec$weights)

# Voltage Fluctuation (Voltage fluctuation days per month)
model.dd.hfe.quality_fluc.scst = fixest::feols(quality_fluc ~ scst*wave + hhsize +
                                                 lnexp + bplaay | hh, 
                                               data = pool_elec, weights = pool_elec$weights)

# Low Voltage (Low Voltage days per month)
model.dd.hfe.quality_low.scst = fixest::feols(quality_low ~ scst*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool_elec, weights = pool_elec$weights)

## Model and coefficient lists for village and household fixed effects combined ##
models.vfe.dd_elec_scst <- list(model.dd.vfe.duration_day.scst,model.dd.vfe.duration_night.scst,
                                model.dd.vfe.reliability.scst, model.dd.vfe.quality_fluc.scst,
                                model.dd.vfe.quality_low.scst)

models.hfe.dd_elec_scst <- list(model.dd.hfe.duration_day.scst, 
                                model.dd.hfe.duration_night.scst, model.dd.hfe.reliability.scst, 
                                model.dd.hfe.quality_fluc.scst, model.dd.hfe.quality_low.scst)

# Generate latex tables
esttex(models.vfe.dd_elec_scst, 
       se = "cluster", 
       digits = 3,
       dict = c(scst = "SC/ST", wave = "2018", duration_day = "Daily Supply Hours",
                duration_night = "Evening Supply Hours", reliability = "Monthly Outage Days",
                quality_fluc = "Fluctuating Voltage Days", quality_low = "Low Voltage Days",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_elec_vfe_scst.tex"))

esttex(models.hfe.dd_elec_scst, 
       se = "cluster", 
       digits = 3,
       dict = c(scst = "SC/ST", wave = "2018", duration_day = "Daily Supply Hours",
                duration_night = "Evening Supply Hours", reliability = "Monthly Outage Days",
                quality_fluc = "Fluctuating Voltage Days", quality_low = "Low Voltage Days",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_elec_hfe_scst.tex"))

# With household fixed effects and two interaction terms -----------------------

# Duration Day (Total hours per day)
model.dd.hfe.duration_day.sc_st = fixest::feols(duration_day ~ sc*wave + st*wave + hhsize +
                                                 lnexp + bplaay | hh, 
                                               data = pool_elec, weights = pool_elec$weights)

# Duration Night (Hours in the evening)
model.dd.hfe.duration_night.sc_st = fixest::feols(duration_night ~ sc*wave + st*wave + hhsize +
                                                   lnexp + bplaay | hh, 
                                                 data = pool_elec, weights = pool_elec$weights)

# Reliability (Outage days per month)
model.dd.hfe.reliability.sc_st = fixest::feols(reliability ~ sc*wave + st*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool_elec, weights = pool_elec$weights)

# Voltage Fluctuation (Voltage fluctuation days per month)
model.dd.hfe.quality_fluc.sc_st = fixest::feols(quality_fluc ~ sc*wave + st*wave + hhsize +
                                                 lnexp + bplaay | hh, 
                                               data = pool_elec, weights = pool_elec$weights)

# Low Voltage (Low Voltage days per month)
model.dd.hfe.quality_low.sc_st = fixest::feols(quality_low ~ sc*wave + st*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool_elec, weights = pool_elec$weights)

models.hfe.dd_elec_sc_st <- list(model.dd.hfe.duration_day.sc_st, 
                                model.dd.hfe.duration_night.sc_st, model.dd.hfe.reliability.sc_st, 
                                model.dd.hfe.quality_fluc.sc_st, model.dd.hfe.quality_low.sc_st)

#Generate latex tables
esttex(models.hfe.dd_elec_sc_st, 
       se = "cluster", 
       digits = 3,
       dict = c(sc = "SC", st = "ST", wave = "2018", duration_day = "Daily Supply Hours",
                duration_night = "Evening Supply Hours", reliability = "Monthly Outage Days",
                quality_fluc = "Fluctuating Voltage Days", quality_low = "Low Voltage Days",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_elec_hfe_sc_st.tex"))

# Visualisation of linear combination of interaction terms ---------------------

# Replicates the stata lincom command, using robust clustered standard errors
# for our specific case of one interaction term (wave) and two binary 
# variables (sc and st) interacted separately

lincom_custom <- function (fit, interactionterm, var1, var2) {
  
  sumvar1interaction <- fit$coefficients[[interactionterm]] + fit$coefficients[[var1]]
  sumvar2interaction <- fit$coefficients[[interactionterm]] + fit$coefficients[[var2]]
  
  vcov <- summary(fit)$cov.scaled
  
  sevar1interaction <- sqrt(
    fit[["se"]][interactionterm]^2 +
      fit[["se"]][var1]^2 +
      2 * vcov[interactionterm, var1])
  
  sevar2interaction <- sqrt(
    fit[["se"]][interactionterm]^2 +
      fit[["se"]][var2]^2 +
      2 * vcov[interactionterm, var2])
  
  lincom_df <- tibble(.rows = 3,
                      "term" = c(interactionterm,var1,var2),
                      "coef" = c(fit$coefficients[[interactionterm]],
                                 sumvar1interaction,sumvar2interaction),
                      "se" = c(fit[["se"]][[interactionterm]],
                               sevar1interaction, sevar2interaction))
  
  return(lincom_df)
  
}

durationday_lincom <- lincom_custom(model.dd.hfe.duration_day.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "Daily Supply Hours",
         type = "Hours")

durationeve_lincom <- lincom_custom(model.dd.hfe.duration_night.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "Evening Supply Hours",
         type = "Hours")

outage_lincom <- lincom_custom(model.dd.hfe.reliability.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "Outage Days",
         type = "Days")

fluctvolt_lincom <- lincom_custom(model.dd.hfe.quality_fluc.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "Fluctuation Voltage Days",
         type = "Days")

lowvolt_lincom <- lincom_custom(model.dd.hfe.quality_low.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "Low Voltage Days",
         type = "Days")

lincom_df <- rbind(rbind(rbind(durationday_lincom, durationeve_lincom),
                         rbind(outage_lincom, fluctvolt_lincom)),
                   lowvolt_lincom) %>% 
  mutate(term = case_when(
    grepl(term, pattern = "sc") ~ "SC",
    grepl(term, pattern = "st") ~ "ST",
    TRUE ~ "All Others"
  ))

lincom_df %>% 
  ggplot(aes(model, coef, colour = term, shape = term,
             ymin = coef - 1.96 * se,
             ymax = coef + 1.96 * se)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(position = position_dodge(width = 0.5), width = 0.5) +
  geom_hline(yintercept = 0, color = "red", linetype = 2) +
  labs(x = "Outcome Variable",
       y = "Absolute change in the outcome variable",
       colour = "Social Group",
       shape = "Social Group"
  ) +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"))

ggsave(here("Manuscript","Figures","improvement_elec_scst.png"),
       device = "png", width = 9, height = 9*0.618, units = "in")
